%随机改变体系构型，选取除原点处的原子外的其他所有原子，产生方向随机，大小为dr * rand的位移
%原子位置限制在r\in[r_min,r_max]的球壳内
%config0是一个1*N的cell（N为体系的原子数），每个元素是一个1*3数组，记录原子坐标
%choice=0时，只改变N>13（外壳层的原子）的位形；choice=1时，改变除原点处原子外的其他所有原子
function config = change_config_multiple(config0,dr, r_max, r_min, choice)
    d = 0.88;
    N = length(config0);
    config  = config0;
%     r_min = 0.95;
%     r_max = 1.15;

    if(choice == 0)
        list = 14:N;
    else
        list = 2:N;
    end
    
    for i = list
        pos_old = config0{i};
        for k=1:1000
            r = dr * rand;
            theta = pi * rand;
            phi = 2*pi * rand;
            dx = r * sin(theta) * cos(phi);
            dy = r * sin(theta) * sin(phi);
            dz = r * cos(theta);
            
            pos_new = pos_old + [dx,dy,dz];
            
            %如果产生的新位置超出空间边界，采用周期边条件的方式处理
            Ri = sqrt(sum(pos_new.^2));
            if(Ri > r_max)
                r_new = Ri - (r_max - r_min);
                x_new = r_new * sin(theta) * cos(phi);
                y_new = r_new * sin(theta) * sin(phi);
                z_new = r_new * cos(theta); 
                pos_new = [x_new, y_new, z_new];
            end
            
            judge = 1;

            for j=1:(i-1)
                rj_vector = config{j} - pos_new;
                rj_square = sum(rj_vector.^2);
                rj = sqrt(rj_square);
                if(rj < d)
                    judge = 0;
                    break
                end
            end

            if(judge == 1)
                config{i} = pos_new;
                break
            else
                config{i} = pos_new;
            end

        end
    end
end



